library(cregg)
library(ggplot2)
library(ggpubr)
rm(list=ls())

#### Load data 

load("data.RData")

#### Formulas for models

f2 <- better_pol ~ age + prof+career + gender + legeff + voteror + voterknow + smedia ## Formula
f3 <- prob_vote ~ age + prof+career + gender + legeff + voteror + voterknow + smedia ## Formula

#### Compute marginal means and AMCEs
marginalmeans.better_pol <- mm(data[data$country!= "Germany" & data$country!="Switzerland",], f2, id = ~id, family = binomial(link = "logit"))
amces.better_pol <- cj(data[data$country!= "Germany" & data$country!="Switzerland",], f2, id = ~id)

marginalmeans.prob_vote <- mm(data[data$country!= "Germany" & data$country!="Switzerland",], f3, id = ~id, family = binomial(link = "logit"))
amces.prob_vote <- cj(data[data$country!= "Germany" & data$country!="Switzerland",], f3, id = ~id)


#### Create data for Figure A1 


dataPlotMarginalmeans.better_pol <- as.data.frame(list(outcome = c(rep("vote_choice", 30)),
                                            feature = c(rep("age", 6), rep("prof", 6), rep("career", 3), rep("gender", 3),
                                                        rep("legeff", 3), rep("voteror", 3), rep("voterknow", 3), rep("smedia", 3)),
                                            level = c("Age", as.character(marginalmeans.better_pol[1:5,]$level),
                                                      "Profession", "Buisness Person", "Farmer", "Lawyer", "Medical doctor", "Parliamentary assistant",
                                                      "Career", "Joined party 3 years ago", "Party offices in last 15 years",
                                                      "Gender", "Man", "Woman",
                                                      "Legislative Effectiveness", "Most projects accepted", "Some projects accepted",
                                                      "Presence in Constituency", "Little with voters in constituency", "Lot of time with voters in constituency",
                                                      "Knowledge of voters", "Low knowledge of voters' position", "Good knowledge of voters' position",
                                                      "Social Media", "Regularly communicate views", "Do not communicate views"),
                                            estimate = c(NA, marginalmeans.better_pol[1:5,]$estimate,
                                                         NA, marginalmeans.better_pol[6:10,]$estimate,
                                                         NA, marginalmeans.better_pol[11:12,]$estimate,
                                                         NA, marginalmeans.better_pol[13:14,]$estimate,
                                                         NA, marginalmeans.better_pol[15:16,]$estimate,
                                                         NA, marginalmeans.better_pol[17:18,]$estimate,
                                                         NA, marginalmeans.better_pol[19:20,]$estimate,
                                                         NA, marginalmeans.better_pol[21:22,]$estimate),
                                            lower = c(NA, marginalmeans.better_pol[1:5,]$lower,
                                                      NA, marginalmeans.better_pol[6:10,]$lower,
                                                      NA, marginalmeans.better_pol[11:12,]$lower,
                                                      NA, marginalmeans.better_pol[13:14,]$lower,
                                                      NA, marginalmeans.better_pol[15:16,]$lower,
                                                      NA, marginalmeans.better_pol[17:18,]$lower,
                                                      NA, marginalmeans.better_pol[19:20,]$lower,
                                                      NA, marginalmeans.better_pol[21:22,]$lower),
                                            upper = c(NA, marginalmeans.better_pol[1:5,]$upper,
                                                      NA, marginalmeans.better_pol[6:10,]$upper,
                                                      NA, marginalmeans.better_pol[11:12,]$upper,
                                                      NA, marginalmeans.better_pol[13:14,]$upper,
                                                      NA, marginalmeans.better_pol[15:16,]$upper,
                                                      NA, marginalmeans.better_pol[17:18,]$upper,
                                                      NA, marginalmeans.better_pol[19:20,]$upper,
                                                      NA, marginalmeans.better_pol[21:22,]$upper)))




dataPlotMarginalmeans.better_pol$level <- factor(dataPlotMarginalmeans.better_pol$level, level = c(dataPlotMarginalmeans.better_pol[19:21,]$level, ## legeff
                                                                             dataPlotMarginalmeans.better_pol[22:24,]$level, ## Voteror
                                                                             dataPlotMarginalmeans.better_pol[25:27,]$level, ## Voterknow
                                                                             dataPlotMarginalmeans.better_pol[28:30,]$level,
                                                                             dataPlotMarginalmeans.better_pol[1:6,]$level, ## Age
                                                                             dataPlotMarginalmeans.better_pol[16:18,]$level, ## Gender
                                                                             dataPlotMarginalmeans.better_pol[13:15,]$level, ## Plo career
                                                                             dataPlotMarginalmeans.better_pol[7:12,]$level)) ## social media

dataPlotMarginalmeans.better_pol$feature <- factor(dataPlotMarginalmeans.better_pol$feature, level = c(unique(dataPlotMarginalmeans.better_pol$feature)[5:8],
                                                                                 unique(dataPlotMarginalmeans.better_pol$feature)[1], 
                                                                                 unique(dataPlotMarginalmeans.better_pol$feature)[4],
                                                                                 unique(dataPlotMarginalmeans.better_pol$feature)[3],
                                                                                 unique(dataPlotMarginalmeans.better_pol$feature)[2]))

dataPlotMarginalmeans.better_pol$type <- "Marginal Means"



dataPlotAmces.better_pol <- as.data.frame(list(outcome = c(rep("vote_choice", 30)),
                                    feature = c(rep("age", 6), rep("prof", 6), rep("career", 3), rep("gender", 3),
                                                rep("legeff", 3), rep("voteror", 3), rep("voterknow", 3), rep("smedia", 3)),
                                    level = c("Age", as.character(amces.better_pol[1:5,]$level),
                                              "Profession", "Buisness Person", "Farmer", "Lawyer", "Medical doctor", "Parliamentary assistant",
                                              "Career", "Joined party 3 years ago", "Party offices in last 15 years",
                                              "Gender", "Man", "Woman",
                                              "Legislative Effectiveness", "Most projects accepted", "Some projects accepted",
                                              "Presence in Constituency", "Little with voters in constituency", "Lot of time with voters in constituency",
                                              "Knowledge of voters", "Low knowledge of voters' position", "Good knowledge of voters' position",
                                              "Social Media", "Regularly communicate views", "Do not communicate views"),
                                    estimate = c(NA, amces.better_pol[1:5,]$estimate,
                                                 NA, amces.better_pol[6:10,]$estimate,
                                                 NA, amces.better_pol[11:12,]$estimate,
                                                 NA, amces.better_pol[13:14,]$estimate,
                                                 NA, amces.better_pol[15:16,]$estimate,
                                                 NA, amces.better_pol[17:18,]$estimate,
                                                 NA, amces.better_pol[19:20,]$estimate,
                                                 NA, amces.better_pol[21:22,]$estimate),
                                    lower = c(NA, amces.better_pol[1:5,]$lower,
                                              NA, amces.better_pol[6:10,]$lower,
                                              NA, amces.better_pol[11:12,]$lower,
                                              NA, amces.better_pol[13:14,]$lower,
                                              NA, amces.better_pol[15:16,]$lower,
                                              NA, amces.better_pol[17:18,]$lower,
                                              NA, amces.better_pol[19:20,]$lower,
                                              NA, amces.better_pol[21:22,]$lower),
                                    upper = c(NA, amces.better_pol[1:5,]$upper,
                                              NA, amces.better_pol[6:10,]$upper,
                                              NA, amces.better_pol[11:12,]$upper,
                                              NA, amces.better_pol[13:14,]$upper,
                                              NA, amces.better_pol[15:16,]$upper,
                                              NA, amces.better_pol[17:18,]$upper,
                                              NA, amces.better_pol[19:20,]$upper,
                                              NA, amces.better_pol[21:22,]$upper)))




dataPlotAmces.better_pol$level <- factor(dataPlotAmces.better_pol$level, level = c(dataPlotAmces.better_pol[19:21,]$level, ## legeff
                                                             dataPlotAmces.better_pol[22:24,]$level, ## Voteror
                                                             dataPlotAmces.better_pol[25:27,]$level, ## Voterknow
                                                             dataPlotAmces.better_pol[28:30,]$level,
                                                             dataPlotAmces.better_pol[1:6,]$level, ## Age
                                                             dataPlotAmces.better_pol[16:18,]$level, ## Gender
                                                             dataPlotAmces.better_pol[13:15,]$level, ## Plo career
                                                             dataPlotAmces.better_pol[7:12,]$level)) ## social media

dataPlotAmces.better_pol$feature <- factor(dataPlotAmces.better_pol$feature, level = c(unique(dataPlotAmces.better_pol$feature)[5:8],
                                                                 unique(dataPlotAmces.better_pol$feature)[1], 
                                                                 unique(dataPlotAmces.better_pol$feature)[4],
                                                                 unique(dataPlotAmces.better_pol$feature)[3],
                                                                 unique(dataPlotAmces.better_pol$feature)[2]))

dataPlotAmces.better_pol$type <- "Average Marginal Component Effect"


png("Figures/Figure A1.png", width = 3000, height = 3000, res=300)
ggplot(rbind(dataPlotAmces.better_pol, dataPlotMarginalmeans.better_pol))+
  geom_point(aes(x=estimate, y=level, color = feature))+
  geom_errorbarh(aes(y=level, xmin=lower, xmax=upper, color=feature), size=0.5, height=0.5)+
  facet_wrap(~type,
             labeller = labeller(
               type = c(`Marginal Means` = "Marginal Means", `amces` = "Average Marginal Component Effect")
             ), scales = 'free_x')+
  ylab("")+
  scale_color_manual(values = c("blue", "green", "black", "grey", "orange", "red", "dark green", "brown"), 
                     labels = c("age" = "Age", "prof"="Profession", 
                                "career"="Career", "gender"="Gender", 
                                "legeff"="Legislative Efficiency", "voteror"="Presence in constituency", 
                                "voterknow"="Knowledge of Voters", "smedia"="Social Media"))+
  geom_vline(aes(xintercept = c(rep(0, 30), rep(.5, 30))), linewidth=.1)+
  theme_minimal()+
  scale_y_discrete(limits = rev)+
  theme(legend.position = "bottom", legend.title = element_blank(), 
        panel.border = element_rect(colour = "black", fill=NA, linewidth=.5), 
        axis.text.y = element_text(face = c(rep("plain", 5), "bold",
                                            rep("plain", 2), "bold",
                                            rep("plain", 2), "bold",
                                            rep("plain", 5), "bold",
                                            rep("plain", 2), "bold",
                                            rep("plain", 2), "bold",
                                            rep("plain", 2), "bold",
                                            rep("plain", 2), "bold")))
dev.off()



#### Create data for Figure A2


dataPlotMarginalmeans.prob_vote <- as.data.frame(list(outcome = c(rep("vote_choice", 30)),
                                                       feature = c(rep("age", 6), rep("prof", 6), rep("career", 3), rep("gender", 3),
                                                                   rep("legeff", 3), rep("voteror", 3), rep("voterknow", 3), rep("smedia", 3)),
                                                       level = c("Age", as.character(marginalmeans.prob_vote[1:5,]$level),
                                                                 "Profession", "Buisness Person", "Farmer", "Lawyer", "Medical doctor", "Parliamentary assistant",
                                                                 "Career", "Joined party 3 years ago", "Party offices in last 15 years",
                                                                 "Gender", "Man", "Woman",
                                                                 "Legislative Effectiveness", "Most projects accepted", "Some projects accepted",
                                                                 "Presence in Constituency", "Little with voters in constituency", "Lot of time with voters in constituency",
                                                                 "Knowledge of voters", "Low knowledge of voters' position", "Good knowledge of voters' position",
                                                                 "Social Media", "Regularly communicate views", "Do not communicate views"),
                                                       estimate = c(NA, marginalmeans.prob_vote[1:5,]$estimate,
                                                                    NA, marginalmeans.prob_vote[6:10,]$estimate,
                                                                    NA, marginalmeans.prob_vote[11:12,]$estimate,
                                                                    NA, marginalmeans.prob_vote[13:14,]$estimate,
                                                                    NA, marginalmeans.prob_vote[15:16,]$estimate,
                                                                    NA, marginalmeans.prob_vote[17:18,]$estimate,
                                                                    NA, marginalmeans.prob_vote[19:20,]$estimate,
                                                                    NA, marginalmeans.prob_vote[21:22,]$estimate),
                                                       lower = c(NA, marginalmeans.prob_vote[1:5,]$lower,
                                                                 NA, marginalmeans.prob_vote[6:10,]$lower,
                                                                 NA, marginalmeans.prob_vote[11:12,]$lower,
                                                                 NA, marginalmeans.prob_vote[13:14,]$lower,
                                                                 NA, marginalmeans.prob_vote[15:16,]$lower,
                                                                 NA, marginalmeans.prob_vote[17:18,]$lower,
                                                                 NA, marginalmeans.prob_vote[19:20,]$lower,
                                                                 NA, marginalmeans.prob_vote[21:22,]$lower),
                                                       upper = c(NA, marginalmeans.prob_vote[1:5,]$upper,
                                                                 NA, marginalmeans.prob_vote[6:10,]$upper,
                                                                 NA, marginalmeans.prob_vote[11:12,]$upper,
                                                                 NA, marginalmeans.prob_vote[13:14,]$upper,
                                                                 NA, marginalmeans.prob_vote[15:16,]$upper,
                                                                 NA, marginalmeans.prob_vote[17:18,]$upper,
                                                                 NA, marginalmeans.prob_vote[19:20,]$upper,
                                                                 NA, marginalmeans.prob_vote[21:22,]$upper)))




dataPlotMarginalmeans.prob_vote$level <- factor(dataPlotMarginalmeans.prob_vote$level, level = c(dataPlotMarginalmeans.prob_vote[19:21,]$level, ## legeff
                                                                                                   dataPlotMarginalmeans.prob_vote[22:24,]$level, ## Voteror
                                                                                                   dataPlotMarginalmeans.prob_vote[25:27,]$level, ## Voterknow
                                                                                                   dataPlotMarginalmeans.prob_vote[28:30,]$level,
                                                                                                   dataPlotMarginalmeans.prob_vote[1:6,]$level, ## Age
                                                                                                   dataPlotMarginalmeans.prob_vote[16:18,]$level, ## Gender
                                                                                                   dataPlotMarginalmeans.prob_vote[13:15,]$level, ## Plo career
                                                                                                   dataPlotMarginalmeans.prob_vote[7:12,]$level)) ## social media

dataPlotMarginalmeans.prob_vote$feature <- factor(dataPlotMarginalmeans.prob_vote$feature, level = c(unique(dataPlotMarginalmeans.prob_vote$feature)[5:8],
                                                                                                       unique(dataPlotMarginalmeans.prob_vote$feature)[1], 
                                                                                                       unique(dataPlotMarginalmeans.prob_vote$feature)[4],
                                                                                                       unique(dataPlotMarginalmeans.prob_vote$feature)[3],
                                                                                                       unique(dataPlotMarginalmeans.prob_vote$feature)[2]))

dataPlotMarginalmeans.prob_vote$type <- "Marginal Means"



dataPlotAmces.prob_vote <- as.data.frame(list(outcome = c(rep("vote_choice", 30)),
                                               feature = c(rep("age", 6), rep("prof", 6), rep("career", 3), rep("gender", 3),
                                                           rep("legeff", 3), rep("voteror", 3), rep("voterknow", 3), rep("smedia", 3)),
                                               level = c("Age", as.character(amces.prob_vote[1:5,]$level),
                                                         "Profession", "Buisness Person", "Farmer", "Lawyer", "Medical doctor", "Parliamentary assistant",
                                                         "Career", "Joined party 3 years ago", "Party offices in last 15 years",
                                                         "Gender", "Man", "Woman",
                                                         "Legislative Effectiveness", "Most projects accepted", "Some projects accepted",
                                                         "Presence in Constituency", "Little with voters in constituency", "Lot of time with voters in constituency",
                                                         "Knowledge of voters", "Low knowledge of voters' position", "Good knowledge of voters' position",
                                                         "Social Media", "Regularly communicate views", "Do not communicate views"),
                                               estimate = c(NA, amces.prob_vote[1:5,]$estimate,
                                                            NA, amces.prob_vote[6:10,]$estimate,
                                                            NA, amces.prob_vote[11:12,]$estimate,
                                                            NA, amces.prob_vote[13:14,]$estimate,
                                                            NA, amces.prob_vote[15:16,]$estimate,
                                                            NA, amces.prob_vote[17:18,]$estimate,
                                                            NA, amces.prob_vote[19:20,]$estimate,
                                                            NA, amces.prob_vote[21:22,]$estimate),
                                               lower = c(NA, amces.prob_vote[1:5,]$lower,
                                                         NA, amces.prob_vote[6:10,]$lower,
                                                         NA, amces.prob_vote[11:12,]$lower,
                                                         NA, amces.prob_vote[13:14,]$lower,
                                                         NA, amces.prob_vote[15:16,]$lower,
                                                         NA, amces.prob_vote[17:18,]$lower,
                                                         NA, amces.prob_vote[19:20,]$lower,
                                                         NA, amces.prob_vote[21:22,]$lower),
                                               upper = c(NA, amces.prob_vote[1:5,]$upper,
                                                         NA, amces.prob_vote[6:10,]$upper,
                                                         NA, amces.prob_vote[11:12,]$upper,
                                                         NA, amces.prob_vote[13:14,]$upper,
                                                         NA, amces.prob_vote[15:16,]$upper,
                                                         NA, amces.prob_vote[17:18,]$upper,
                                                         NA, amces.prob_vote[19:20,]$upper,
                                                         NA, amces.prob_vote[21:22,]$upper)))




dataPlotAmces.prob_vote$level <- factor(dataPlotAmces.prob_vote$level, level = c(dataPlotAmces.prob_vote[19:21,]$level, ## legeff
                                                                                   dataPlotAmces.prob_vote[22:24,]$level, ## Voteror
                                                                                   dataPlotAmces.prob_vote[25:27,]$level, ## Voterknow
                                                                                   dataPlotAmces.prob_vote[28:30,]$level,
                                                                                   dataPlotAmces.prob_vote[1:6,]$level, ## Age
                                                                                   dataPlotAmces.prob_vote[16:18,]$level, ## Gender
                                                                                   dataPlotAmces.prob_vote[13:15,]$level, ## Plo career
                                                                                   dataPlotAmces.prob_vote[7:12,]$level)) ## social media

dataPlotAmces.prob_vote$feature <- factor(dataPlotAmces.prob_vote$feature, level = c(unique(dataPlotAmces.prob_vote$feature)[5:8],
                                                                                       unique(dataPlotAmces.prob_vote$feature)[1], 
                                                                                       unique(dataPlotAmces.prob_vote$feature)[4],
                                                                                       unique(dataPlotAmces.prob_vote$feature)[3],
                                                                                       unique(dataPlotAmces.prob_vote$feature)[2]))

dataPlotAmces.prob_vote$type <- "Average Marginal Component Effect"


png("Figures/Figure A2.png", width = 3000, height = 3000, res=300)
ggplot(rbind(dataPlotAmces.prob_vote, dataPlotMarginalmeans.prob_vote))+
  geom_point(aes(x=estimate, y=level, color = feature))+
  geom_errorbarh(aes(y=level, xmin=lower, xmax=upper, color=feature), size=0.5, height=0.5)+
  facet_wrap(~type,
             labeller = labeller(
               type = c(`Marginal Means` = "Marginal Means", `amces` = "Average Marginal Component Effect")
             ), scales = 'free_x')+
  ylab("")+
  scale_color_manual(values = c("blue", "green", "black", "grey", "orange", "red", "dark green", "brown"), 
                     labels = c("age" = "Age", "prof"="Profession", 
                                "career"="Career", "gender"="Gender", 
                                "legeff"="Legislative Efficiency", "voteror"="Presence in constituency", 
                                "voterknow"="Knowledge of Voters", "smedia"="Social Media"))+
  geom_vline(aes(xintercept = c(rep(0, 30), rep(5, 30))), linewidth=.1)+
  theme_minimal()+
  scale_y_discrete(limits = rev)+
  theme(legend.position = "bottom", legend.title = element_blank(), 
        panel.border = element_rect(colour = "black", fill=NA, linewidth=.5), 
        axis.text.y = element_text(face = c(rep("plain", 5), "bold",
                                            rep("plain", 2), "bold",
                                            rep("plain", 2), "bold",
                                            rep("plain", 5), "bold",
                                            rep("plain", 2), "bold",
                                            rep("plain", 2), "bold",
                                            rep("plain", 2), "bold",
                                            rep("plain", 2), "bold")))
dev.off()



